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Abstract.  In  solving  stiff  systems  of  oidinaiy  differential  equations  using  BDF 
methods,  Jacobians  needed  for  quasi-Newton  itentioo  are  frequently  computed  using 
finite  differences.  Round-<rff  errors  in  the  finite*difference  approximadcm  can  lead  to 
Newton  failures  forcing  the  code  to  choose  its  time  steps  based  on  ^‘stability”  rather 
dian  accuracy  considenuians.  When  standard  stqtsire  control  is  used  the  code  can 
etqierienoe  thrashing  which  increases  the  total  number  of  time  steps,  Jacobian  evalua¬ 
tions,  and  function  evaluations.  In  this  paper  we  investigate  this  situation,  explaining 
some  surprising  time  step  selection  behavior  produced  by  the  standard  control  mechan¬ 
ism.  A  new  control  mechanism  is  proposed  which  attempts  to  find  and  use  a  “stabil- 
ity”  stqnize.  A  comparison  of  the  new  strategy  with  the  standard  strategy  and  with 
two  PI  controllers  introduced  earlier  is  made  using  the  stiff  test  set. 

L  latrodnctioo. 

In  solving  stiff  systems  of  ordinary  differential  equations  using  BDF  methods, 
Jacobians  which  are  needed  for  the  quasi-Newton  iteration  are  ofren  iqiproximated 

*  The  worit  of  tbit  tabu  was  paniaUy  lupponed  by  ARO  CootEKt  Nonaber  DAAL03-89-C-0038 
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using  finitt  diHerences  [1,  2,  9).  SmaUer  stepsiaes  than  allowed  by  accuracy  con- 
sideiaiiont  may  be  needed  to  guataniee  convergence  of  the  Newton  iteration  due  to 
round-off  errors  in  the  fihite-difieieace  Jacobians.  The  standard  stepsize  control 
mechanism,  such  as  that  used  in  DASSL  [2],  is 

h,  *  (TOL/£Sr,.i)'^  (1.1) 

udwre  TOL  is  the  user-prescribed  tokranoe,  EST^.i  is  the  local  error  esdmaie  com¬ 
puted  at  time  ^  p  is  the  method  order,  fhnvever,  (1.1)  is  based  solely  on  accu¬ 
racy  considerations.  This  can  lead  to  hi^y  oscillatory  stepsize  behavior  (see  Figure 
1.1). 

Here  we  apply  DASSL  to  the  stiff  test  set  [4,  S]  with  an  approximate  Jact^ian  to 
simulate  the  effects  of  an  inaccurate  finite  difference  Jacobian.  The  stepsize  behavior 
shown  in  Figure  1.1  when  DASSL  is  applied  to  problem  D2  is  typical.  After  periods 
of  taking  relatively  small  stepsizes  dx;  algorithm  suddenly  increases  the  stepsize  by 
several  orders  of  magnitude.  It  remains  at  this  larger  stepsize  for  several  time  steps, 
and  then  decreases  the  stepsize  dramatically  whereupon  the  process  begins  again. 
Although  most  of  the  tune  steps  taken  by  DASSL  are  of  the  smaller  size,  the  solutitm 
on  most  cf  the  time  interval  is  found  using  the  larger  time  steps. 

In  92  we  analyze  the  behavior  of  DASSL  oa  a  simple  linear  system  which  leads 
to  an  understanding  of  the  stepsize  briutvior  discusated  above.  We  also  present  a 
mndificatinn  of  the  time  step  selection  strategy  used  in  DASSL  based  on  the  quasi- 
Newion  algorithm  of  Dennis  and  Schnabel  [3].  The  revised  strategy  prevents  the 
larger  “anomalous"  steps  but  leads  to  a  much  larger  number  of  rime  steps  with  no 
significant  in^uoyement  in  accuracy.  The  stepsize  ccMitroiler  of  Gustafsson  et  al.  [6], 
referred  to  hencdbcth  as  the  PI  cootrcriler  L  is  presented  in  93  with  a  new  interpreta¬ 
tion.  This  contrtriler  was  developed  for  explicit  time  integrators.  Gustafeson  [8] 
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(teveloped  a  n  controller  (tefetred  to  lurein  as  PI  ctmtroUer  II)  for  implicit  methods. 
We  also  discuss  this  controller  in  §3.  In  §4  we  present  a  new  control  strategy  and 
cooqMue  it  with  the  standard  stepsize  controller  and  the  PI  controllers  on  the  stiff  test 
set  (4, 5].  Brief  ccmclusioos  are  given  in  §5. 


2.  Anniysb  of  a  stiff  system. 

Consider  ISrst  the  standard  test  problem 

}''  =  Xy.  y(0)=l  (2.1) 

where  Re  (k)  ^  0.  Applying  the  backward-Euler  method  together  with  quasi-Newton 
iteration  yields 

(1  -  aXA)Ay‘^,  =  -  XA)  +y,.  i  ^  1.  (2.2) 

Typically,  if  an  analytic  Jacobian  is  used  a  s  1.  To  model  the  effects  of  an  inaccurate 
matrix  approximation,  we  choose  a  different  fn»n  1.  After  k+l  iterations  of  the 
quasi-Newton  method  we  obtain 

_  0  y.  (g-DXA, 

yn*i  (i_aXA)*-^‘  l-oXA^  l-oXA 


(a-l)*(-XA)* 

(l-oXA)* 


(2.3) 


where  yJVi  is  the  predicted  solution.  If  iXA(a  -  1)/(1  -  aXA)i  <  1  the  quasi-Newuxi 
iterates  converge  to  the  the  true  solution  where  the  rate  of  convergence,  p.  is  given  by 

-y«/(l  -XA)I 


^  =  lXA(a-  1)/(1  -aXA)l. 

-va-Wi)l 

In  DASSL  [2]  the  quasi-Newton  iteration  is  said  to  converge  if 


(2.4) 


l-() 


W;'  -J’f.il<0.33 


(2j») 
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wliere  0  is  an  apfsoximatioo  ci  the  rate  p  given  by 

-  (lyitl  -  yt-.,  i'iyn*i  -  y.i,  i)".  (2JW 

Thus  die  numbo*  of  iterations  before  convergence  is  determined  by  both  the  accuracy 

the  predictor  and  the  rate  of  convergence  0. 

For  systems  of  equations,  the  analysis  is  complicated  by  the  nom  used. 
Aldiough  DASSL  uses  a  weighted  rms  norm,  herein  we  consider  the  norm  which 
displays  the  same  type  of  bdiavior.  Consider  the  diagonal  system 

/*Dy.  y(0)=l,  (2.6) 

where  D  *  diag(X],Xo,  X^)  with  <  0,  /  =  1,  2,  •••,  m.  After  two  iterations 

the  rate  of  cmivergence  is  given  by 

[(PiO'.Vi  -y..i/a 

I(y “.1,1  -  y..i'0  -  M))“ + -  +  Ohinj.  - 

where  tire  second  subscript  indicates  the  component  of  the  vector  y  and  p,-  is  given  by 
(2.4)  with  the  apprt^jtiate  X,-.  Unlike  (2.1)  the  rate  of  convergence  of  (2.7)  is  not  con¬ 
stant,  but  is  instead  determiired  by  the  stiffness  (through  p,- )  and  the  accuracy  of  the 
pi'Bflictor  -  Xf/i)).  Thus  if  the  stiff  components  axe  sufficiently  mote 

accurately  predicted  than  the  nonsdff  congxinents  (which  is  likely  since  the  stiff  com¬ 
ponents  change  litde  from  stqi-io-step)  the  rate  of  convergence  will  be  controlled  by 
the  rate  of  convergence  for  die  nonsdff  conqxnients  which  have  snudler  rate  constants. 
The  stepsize  controller  may  then  wish  to  increase  the  stepsize  (which  may  be  low  for 
the  mmstiff  components)  until  the  errors  in  the  stiff  components  are  excited  whereupon 
the  stepsize  undergoes  a  drastic  reduction  since  the  rate  is  now  being  determined  by 
die  stiff  components. 

That  this  actually  luq^rens  can  be  seen  by  applying  DASSL  to  problem  A4  of  the 
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s«lff  test  set  [5]  which  has  ki  *  i  «  I,  2,  •••,  10.  We  soived  this  probiem  on 
0  <  t  0.4  with  absolute  and  relative  error  tolerances  of  0.01  and  with  a  =  0.3.  Fig¬ 
ure  2.1  shows  the  time  steps  and  the  error  in  the  stiffest  component  over  the  interval. 
When  the  stiff  cmxiponent  becomes  sufficiently  accurate  (after  very  r,maU  time  steps) 
the  time  step  increases  rapidly,  reaching  a  value  controlled  by  the  error  in  nonstiff 
cocnpooents.  Such  large  time  steps  excite  errors  in  the  stiff  components  until  the  rate 
of  convergence  becomes  dominated  by  the  stiff  comptments  and  the  nme  step  is  drasti¬ 
cally  reduced,  beginning  the  process  again.  Such  behavior  is  also  seen  in  non-diagonal 
systems  as  shown  in  Figure  1.1. 

The  quesdon  arises,  is  it  desirable  to  permit  the  large  stepsizes.  We  observed  that 
on  some  successful  steps  at  the  larger  stepsizes  IlgCr.Y.YOH  increased  where  Y  is  the 
BDF  solution  when  we  are  solving 

g(tjy)«0.  (2.8) 

One  approach  we  implemented  to  omtect  this  problem  was  to  insist  that  iig(r,Y,Y')]i 
decrease  before  a  step  was  converged.  Specifically  we  required  that  on  each  Newton 
step 

J|gYl.4y)lP-|lt(Y)f  2  0.05  (2.9) 

-2|l«(Y)|f 

(cf.  Dennis  and  Schnabel  [3])  where  AY  is  the  quasi-Newton  direction.  Although 
using  (2.9)  did  reduce  the  number  and  extent  of  the  large  stepsize  regions,  more  time 
steps  were  used  with  no  significant  improvement  in  accuracy.  Thus,  it  seems  reason¬ 
able  to  allow  the  large  steps  and  subsequently  we  do  not  use  (2.9). 

3.  PlcoiiCroL 

As  wu  seen  in  §2  the  standard  stepsize  ctuitrol  mechanism  (1.1)  leads  to  osciila- 
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tcvy  time  ste{».  The  difficulty  is  finding  a  way  to  smooth  out  the  selection  of  small 
time  steps  without  eliminaring  the  selection  of  the  large  time  steps.  One  possible 
approach  is  to  use  PI  controller  I  introduced  by  Gustafsson  et  aL,  [6,  7,  8]  for  explicit 
methods  or  PI  ctmtroller  II  of  Gustafsson  [8]  for  implicit  methods.  In  this  section  we 
present  a  new  interpretation  of  these  strategks  along  with  some  modifications  for  use 
in  our  situation. 

PI  controller  I  can  be  written  in  die  form  (with  some  modifications  for  maximum 
rate  of  stepsize  increase  and  decrease)  [7] 

/i„  =  (TOLlEST^.0‘^(ESr^.2/EST^-f^h^.i  (3.1) 

where  £ST„.2  and  £57^.)  are  estimates  of  the  local  truncation  error  at  time  ^ 
r,«i,  respecdvely.  and  K/  and  Kp  are  parameters  whose  values  depend  only  on 
whedier  a  step  is  successful  or  not  Values  for  AT/  and  Kp  are  given  in  [6]  and  [7] 
aldioagh  diey  differ  slightly.  For  our  purposes  it  is  important  to  note  first  that  in  the 
case  of  a  rejected  step  AT/  Up  and  Kp  -  0.  Thus,  when  a  step  is  rejected  the  stan¬ 
dard  controller  is  used.  Second,  in  the  case  of  an  accepted  step  AT/  +  Kp  =■  Mp . 

We  can  rewrite  (3.1)  as 

K  *  (£Sr,.2/7’OL)^'’(rGL/£Sr„_i)^^'  "  (3.2) 

Assuming  Ki  +  Kp  -  Mp  and  using 

^aeeai-l  (3.3a) 

we  obtain 

K  =  {EST^.2/T0L  (3.3b) 

where  Aoccji-i  represents  the  stepsize  based  on  the  local  truncaticm  error  that  could 
have  been  taken  at  time  r„_2‘  ^  represents  the  stepsize  based 
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00  aocuncy  that  is  typically  used  for  the  next  time  step  Now,  since  (3.3a) 

holds  at  time  1^-2  we  obtain 

K  =  (3.4) 

where  »  KpliKp  +  AT/).  Thus  the  new  sepsize  is  chosen  to  be  the  stepsize  based 
oo  accuracy  multiplied  by  a  factor  that  rqnesents  the  ratio  of  the  actual  stepsize  to 
accuracy  stepsize  we  could  have  chosen  (»  the  prevkxis  stq>.  If  on  the  last  step  the 
accuracy  and  actual  time  step  were  the  same,  the  accuracy  stepsize  is  used  on  the 
present  step.  Otherwise  (the  previous  stepsize  can  never  be  larger  than  the  previous 
stepsize  based  on  accuracy),  if  the  previous  stepsize  was  much  smaller  than  the  previ¬ 
ous  stepsize  based  on  accuracy  only  a  firacdon  of  the  accuracy  stepsize  is  used. 

A  similar  analysis  shows  that  PI  ctxittoller  II  has  the  form 

(3*5) 

if  two  or  more  successive  accepted  time  steps  have  been  taken  (otherwise  the  standard 
controller  is  used).  Now  the  accuracy  time  step  haa;jt.\  is  multiplied  by  an  additional 
factor  representing  the  ratio  of  successive  accepted  time  steps. 

The  version  of  the  PI  algtmthms  we  used  in  our  testing  consists  of  three  cases. 

1:  If  the  present  step  is  rejected  due  to  the  error  test,  set  /i„  =  ^acc,n-i> 

2:  If  the  present  step  is  rejected  due  to  Newton  divergence,  set 

3:  If  the  present  step  is  accepted  update  Kq  s  maxffac^AT^;  JCyj )  if  the  previous  step 
was  not  a  Newton  failure.  For  PI  controller  I  set 
hn  *  min(l,(A„_2//»acc^_2)*")haee.ii-i  set 
hff  *  min(l,(A„_2//la£cji-2) 


for  PI  cmitroller  n. 
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where  fac  *•  0.9,  Kijq  »  OJ,  and  ^Ht  a  0.7  are  fixed  panuneters.  This  differs  slightly 
fiom  the  ^^Koach  taken-  by  Gtistafsson  [6.  8],  since  we  allow  Kq  to  vary.  We  found 
varying  Kq  resulted  in  sli^tly  better  performance  over  fixed  Kq  . 

Gustafsson  [8]  offers  some  addidcnud  ideas  for  PI  controller  11.  He  has  a  stepsize 
algorithm  for  the  case  of  successive  stepsize  failure  due  to  error  control.  In  our  situa¬ 
tion,  however,  we  encounter  successive  stqmze  failure  due  to  divergence  in  Newton’s 
method  so  we  did  not  incorparate  this  heuristic  in  our  algorithm.  In  certain  cases  of 
Newtcm  divergence  he  computes  a  second,  “stability”  stepsize  which  is  based  on  the 
size  of  the  norm  of  the  Jacobian.  Since  we  are  also  interested  in  solving  differential- 
algebraic  equations,  we  are  very  reluctant  to  use  a  scale-dependent  quantity  in  our 
algorithm  so  we  have  neglected  this  feature. 

4.  A  new  controller. 

Fot  reasons  that  will  become  clear  fimn  the  examples  in  this  section  neither  the 
standard  controller  nor  the  PI  controllers  possess  the  desired  stepsize  behavior.  Using 
the  analysis  from  §2  we  present  a  new  controller  which  we  refer  to  as  the  STAB  con¬ 
troller.  We  then  present  a  numerical  comparison  of  the  the  standard,  PI,  and  STAB 
controllers  q;>plied  to  the  stiff  test  set  [4, 5]. 

Bom  our  observations  in  §2  we  desire  a  controller  that  interferes  with  the  stan¬ 
dard  controller  as  little  as  possible.  Our  gmd  was  to  smooth  out  the  time  step  selec¬ 
tion  strategy  <mly  when  the  code  is  thrashing  due  to  Newton  convergence  difficulties. 
When  the  code  is  able  to  use  larger  stepsizes  because  of  a  good  predictor  for  the  stiff 
components,  we  want  to  let  it  do  this  because  this  is  where  it  makes  most  of  its  pro¬ 
gress.  Addititmally,  as  indicated  in  S2,  no  accuracy  is  lost  in  accepting  the  large  steps. 
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We  begin  with  the  observation  that  there  axe  two  important  time  step  sizes,  an 
accuracy  size  and  a  stability  size  where  here,  stability  refers  to  the  conver¬ 
gence  of  Newton’s  method.  Normally  is  smaller  than  /t^  for  BDF  methods  but 
when  the  Jacobian  is  poorly  approximated  the  reverse  can  occur.  In  the  graph  on  the 
right  side  of  Figure  1.1,  the  peaks  in  curve  A  indicate  stepsizes  chosen  by  accuracy  but 
which  caused  the  quasi-Newton  iteration  to  diverge.  Thus  the  best  guess  at  is 
represented  by  the  peaks,  although  it  may  be  quite  a  bit  larger  (curve  C  on  the  right, 
Hgure  1.1).  After  each  peak  two  successful,  smaller  steps  are  taken.  The  value 
is  approximated  by  curve  B  in  Figure  1.1.  Our  controller  seeks  to  detect  when  is 
smaller  than  and  then  makes  two  attempts  at  finding  /t^.  The  algorithm  then 
limits  the  time  step  size  for  10  steps  to  0.87^^^^  after  which  it  reverts  to  the  standard 
ccmtroiler.  The  value  0.87  is  chosen  to  reduce  the  number  of  step  failures  and  to 
reduce  the  number  of  Newton  iterates  required  for  convergence.  If  the  time  step  used 
was  Newton  make  take  several  su^ps  to  converge  due  to  a  larger  rate  and  a 
poorer  predictor. 

The  new  controller  is  invoked  only  when  the  Newton  iteration  fails  to  converge 
(as  long  as  the  convergence  is  not  due  to  a  singular  Jacobian),  i.e.,  when  the  criteria 
(2Ja)  fails  and  when  the  last  successful  time  step  is  smaller  than  the  first  failed 
(Newton)  step  />„.  Two  attempts  are  made  at  finding  hj^.  After  the  first  Newton 
failure,  if  ^  0.8,  h„  =  0.87h„_i  and  the  time  step  is  not  allowed  to  become 

larger  than  this  value  for  10  time  steps  (of  course  it  can  become  smaller  due  to  subse¬ 
quent  Newton  failures  or  error  failure).  If  hn-xlh^  <  0.8  then  h„  =  0.8h„_i  +  0.2/i„. 
Now,  howevo*,  the  stepsize  is  allowed  to  increase,  but  at  a  reduced  maximum  rate  of 
1.18.  If  Newton  fails  for  a  second  (but  not  secmid  consecutive)  time  within  the  10 
step  limit,  /t„  =  0.87hn_i  and  no  increase  above  this  value  is  allowed  for  10  steps. 
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Consecudve  Newton  failures  result  in  tl^  algorithm  reverting  to  the  standard  controller 
since  we  no  longer  seem  to  have  a  good  ^>proximation  to  hg^.  We  limit  our  con¬ 
troller  to  10  steps  so  that  largo-  sta^ize  increases  are  allowed  from  time  to  time  which 
should  preserve  the  desirable  property  of  the  standard  controller. 

We  solved  the  stiff  set  using  DASSL  with  the  three  controllers  and  absolute  and 
relative  error  tolerances  of  10~*,  k  6  and  a  =  0.5.  Table  4.1  contains  the 

number  of  time  steps  used  (including  successful  and  unsuccessful  time  steps)  by  each 
of  the  three  algorithms  with  tolerances  of  0.0001.  None  of  the  algorithms  was  able  to 
solve  FI  or  F4  in  10(X)0  time  steps  and  only  PI  controller  I  was  able  to  do  so  for  FS 
with  this  poor  approximation  to  the  iteration  matrix.  In  alnoost  all  cases  the  STAB 
controller  outperforms  the  standard  controller.  For  problems  A2,  D3,  and  DS  where 
the  standard  controller  appears  to  ou^erform  the  STAB  controller,  the  error  produced 
using  the  STAB  controller  was  at  least  a  factor  of  three  smaller.  Also  note  the 
anomalous  behavior  of  the  standard  ccmtroller  on  F3.  PI  connoUer  11  is  slightly  better 
than  PI  controller  1.  The  STAB  controller  is  mme  efficient  than  the  PI  controllers  for 
a  significant  number  of  the  test  problems.  Although  PI  controller  II  appears  to  per¬ 
form  better  than  the  STAB  controller  cm  the  early  problems  A2,  A3,  and  A4,  the 
STAB  controller  produced  solutions  with  errcns  at  least  a  factor  of  10  smaller  than  PI 
controller  n.  (Dn  problems  D4  and  DS  the  STAB  controller  is  more  accurate  by  a  fac- 
mr  of  at  least  2.5.  The  results  at  this  tolerance  were  typical  of  the  performance  of  the 
ctmtrol  algmithms  at  the  other  tolerances. 

In  Table  4.2  we  have  listed  the  number  of  Jacobian  evaluations,  function  evalua¬ 
tions,  and  errors  (in  the  norm)  for  several  cases.  The  STAB  controller  is  clearly 
more  efficient  in  almost  every  case.  For  several  problems  it  uses  a  factor  of  6  fewer 
Jacobian  evaluations  and  less  than  half  the  number  of  function  evaluations  while 
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obtaining  a  smaller  error.  Even  in  the  cases  where  it  uses  mOTe  function  evaluations 
and  time  steps,  such  as  D4,  its  solution  is  more  accurate.  PI  controller  II  is  generally 
superior  to  PI  controller  I.  The  PI  controllers  have,  in  general,  larger  eirors  and  are 
often  between  the  STAB  and  standard  controllers  in  the  other  measures. 

Hgure  4.1  presents  the  time  steps  used  in  the  standard  controller  and  the  STAB 
controller  for  0  <  r  ^  8  with  tolerances  of  0.00001  and  a  =  0.5.  Qearly  the  STAB 
controller  produces  a  smoother  time  step  history  than  the  standard  controller. 

The  only  problem  which  presented  any  difficulty  to  DASSL  with  the  standard 
controller  and  a  =  1  was  FS.  We  solved  F5  with  the  standard,  the  STAB,  and  the  PI 
ctmtrollers  with  a  =  1  and  with  both  analytic  and  finite  difference  Jacobians  for  toler¬ 
ances  of  10~*,  k  s  2,  3,  6.  The  results  shown  in  Table  4.3a  (analytic  Jacobians) 

and  Table  4.3b  (finite-difference  Jacobians)  indicate  that  the  STAB  controller  enhances 
the  perfoimance  of  DASSL  in  both  cases. 

5.  CtNidicnons. 

When  finite-difference  approximations  to  the  Jacobian  are  used  in  stiff  solvers 
such  as  DASSL,  thrashing  of  the  time  step  can  occur  because  the  accuracy  stepsize 
exceeds  the  stepsize  required  for  conmgence  of  the  quasi-Newton  iteradon.  After 
damping  the  stiff  components  at  small  stepsizes,  such  algorithms  using  the  standard 
stepsize  control  mechanism  are  able  to  take  larger  time  steps  based  on  the  error  in  the 
nonstiff  components.  However,  the  stiff  components  then  become  excited  resulting  in 
drastic  decreases  in  the  stepsize.  We  analyzed  a  simple  linear  system  to  explain  this 
stepsize  phenomena.  To_smooth  out  the  smaller  stepsizes  we  tried  modified  versions 
of  two  PI  controllers  proposed  by  Gustafsson  et  al.  [6,  7,  8]  for  explicit  Runge-Kutta 
methods  and  implicit  Runge-Kutta  methods,  respectively,  which  we  reinterpreted  to  our 
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nouuion.  A  new  controller  was  also  pnqx>sed  which  attempts  to  find  and  use  a  “sta¬ 
bility”  stepsize.  In  comparing  the  four  controllers  on  the  stiff  test  set  with  altered 
Jacotnan  we  found  the  new  controller  superior  in  almost  every  case.  In  fact  for  a 
number  of  cases  it  offers  dramatic  improvement  over  the  standard  and  the  PI  controll¬ 
ers. 


Although  most  of  the  test  results  of  §4  were  somewhat  artificial  the  new  con¬ 
troller  does  produce  significantly  better  results  especially  in  terms  of  the  number  of 
Jacobian  evaluations.  This  continued  to  be  true  even  when  analytic  and  finite 
difference  Jacobians  were  used  for  problem  F5.  Smoothing  out  the  smaller  stepsizes 
also  seems  to  be  beneficial  for  the  error.  More  testing  needs  to  be  done  to  verify  the 
usefulness  of  the  algorithm  in  a  wider  setting,  when  us^  to  solve  partial  differential 
equatimis  by  the  method-of-lines  and  differential-algebraic  equations. 
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Tim  Tim 

Figure  1.1.  Hine  stq)s  used  in  sotving  pwiMem  D2  firooi  the  Hiff  left  set  (5) 
on  0  <  /  ^  8  (left)  and  on  1.83  ^  ^  1.98  (right  dcnoied  by  A)  uaag 
DASSL  with  absolute  and  relative  enor  tolenuioes  of  0.01  and  a  «  OJ.  Es- 
rimate*  of  the  stepsize  based  on  ‘^stabili^”,  (ri^t,  denoted  by  B)  and 
accuracy,  hg^c  (right,  denoted  by  C). 
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Hgnre  2.1.  Time  stqis  (A)  ad  orar  in  7io  (B)  in  loiving  A4  finom  the 
st^  test  set  [S]  on  0  <  f  0.4  widt  DASSL  with  absolute  and  relative 
eiTor  tolerances  of  0.01  and  a  s  0.5. 
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Stiff  Set 

Frobten 

Number  of  Time  Steps 

Standard 

Controller 

STAB 

Controller 

PI 

Controller  I 

PI 

CotHP**!^**  n 

A1 

367 

217 

286 

343 

A2 

1546 

1627 

1990 

984 

A3 

2446 

1467 

1946 

746 

A4 

1469 

1417 

1201 

B1 

1154 

892 

2578 

2372 

B2 

no 

106 

120 

106 

B3 

145 

120 

132 

127 

B4 

277 

184 

360 

328 

BS 

1633 

1071 

1413 

1334 

Cl 

347 

238 

153 

C2 

328 

242 

258 

274 

C3 

322 

152 

257 

303 

'  C4 

348 

281 

420 

383 

C5 

272 

233 

393 

358 

D1 

2099 

621 

1928 

2023 

D2 

2767 

1576 

1810 

D3 

376 

449 

358 

364 

D4 

154 

497 

99 

106 

D5 

144 

182 

131 

107 

D6 

517 

476 

510 

El 

62 

73 

73 

E2 

156 

156 

871 

768 

E3 

2121 

1063 

1973 

1635 

E4 

1298 

972 

1339 

1270 

E5 

.  19 

19 

19 

19 

F2 

159 

135 

112 

138 

F3 

15002 

152 

15002 

15002 

Table  4.1.  The  number  of  time  steps  needed  by  DASSL  with  the  stan- 
daid  controller,  the  STAB  controller,  and  the  PI  conooUers  in  striving 
die  problems  of  the  stiff  test  s^  [4,  5]  with  tolerances  0.0001  and 
a  a  OJ.  None  of  the  three  algoriduns  solved  FI  or  F4  in  fewer  than 
10000  steps  and  only  PI  controller  I  was  able  to  strive  F5  in  fewer  than 
10000  steps. 
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Stiff  Set 
Problem 

TOL 

Standard 

Controller 

STAB 

CjontroUer 

PI 

Controller  I 

PI 

controller  II 

A2 

i(r* 

3349 

8649 

3595 

0.42 

697 

6423 

2632 

0.11x10-6 

2361 

7182 

2918 

0.68x10"^ 

1030 

3632 

1483 

O.MxKT* 

JACS 

FNS 

STEPS 

/2  error 

lOr^ 

8241 

22570 

8630 

0.19  xlOr^ 

1077 

10237 

4110 

0.54  xlOr^ 

2679 

10037 

3958 

0.25  xlOr* 

1117 

6218 

2742 

0.46  xKT^ 

JACS 

FNS 

STEPS 

/2  error 

A3 

ICr^ 

6002 

15458 

4375 

0.36  xir^ 

905 

8554 

2933 

0.43  xlOr^ 

3666 

11065 

4312 

0.35x10-2 

492 

3079 

1375 

0.62x10-2 

JACS 

FNS 
STEPS 
/2  ERROR 

10-^ 

9901 
25928 
10595 
0.96  xKT^ 

1207 

11786 

4594 

0.14xl(r3 

1906 

10087 

4109 

0.20xl(r2 

2883 

10608 

4384 

0.21  X 10-2 

JACS 

FNS 
STEPS 
/2  ERROR 

D2 

i(r^ 

5223 

13226 

5377 

0.11  X 10"^ 

704 

6431 

2607 

0.12xir* 

l092 

10014 

3788 

0.12xl(r‘ 

2988 

8863 

3420 

0.20  xir' 

JACS 

FNS 
STEPS 
/2  ERROR 

i(r* 

9123 

24408 

0.28  xir^ 

967 

9806 

3761 

0.29x10-^ 

4519 

18768 

6762 

0.44xl(r2 

7224 

20259 

7554 

0.51  X  ir2 

JACS 

FNS 

STEPS 

(2  error 

D4 

icr^ 

488 

1129 

506 

0,19  xir^ 

434 

3520 

1538 

0.12xl(r3 

1029 

2500 

1066 

0.64  xltr^ 

311 

824 

341 

0.10x10-2 

JACS 

FNS 
STEPS 
/2  ERROR 

i(r® 

1998 

4859 

2129 

0.18  xlOr^ 

589 

5353 

2330 

0,35x10-* 

209 

1058 

482 

0.30  xlOr^ 

166 

856 

389 

0.30x10-^ 

JACS 

FNS 
STEPS 
/2  ERROR 

E3 

10“^ 

5403 

14102 

5609 

0,21x10-* 

521 

5905 

2120 

0,28x10^ 

3393 

12103 

4293 

0.31  x  lOr^ 

4599 

12439 

4873 

0.23x10-2 

JACS 

FNS 
STEPS 
/2  ERROR 

Table  42.  The  number  of  Jacobian  evaluatkHis  (JACS),  function  evaluations 
(FNS),  time  sfq)s  (STEPS),  and  enor  in  the  norm  for  selected  problems 
foom  the  stiff  test  set  [5]  using  DASSL  with  standard,  STAB,  and  PI  con- 
trdleis. 
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Analytic 


TOL 

Jacobian 

Standard 

Omtioller 

STAB 

Cootndler 

PI 

Controller  I 

PI 

Contiolter  n 

msm 

49 

28 

51 

49 

JACS 

82 

83 

86 

82 

FNS 

49 

44 

52 

49 

39 

51 

84 

39 

JACS 

lOr^ 

71 

170 

142 

71 

FNS 

39 

87 

85 

39 

STEPS 

356 

100 

313 

421 

JACS 

Kr-* 

637 

314 

545 

687 

FNS 

370 

145 

329 

430 

STEPS 

2575 

2137 

JACS 

i(r^ 

4138 

3495 

FNS 

2604 

195 

STEPS 

9296 

3547 

12576 

12211 

JACS 

lor^ 

14798 

12166 

19950 

FNS 

9339 

5675 

12710 

12326 

STEPS 

Table  4.3a.  The  number  of  Jacobian  evalnatimis  (JACS),  function  evalua¬ 
tions  (FNS),  and  dme  steps  (STEPS)  needed  to  solve  problem  F5  firom  the 
stiff  test  set  with  DASSL  and  analytic  Jacobians  using  the  standard,  STAB, 
and  PI  controllers. 


TOL 

Hniie-Diffeience 

Jacobian 

:  Standard 

1  Controller 

STAB 

Controller 

PI 

Controller  I 

PI 

Controller  n 

Tiai 

34 

43 

35 

JACS 

■■ 

226 

248 

202 

FNS 

K9 

47 

44 

35 

STEPS 

74 

48 

39 

74 

JACS 

418 

346 

230 

418 

FNS 

74 

74 

40 

74 

STEPS 

411 

146 

38 

253 

JACS 

lor* 

2328 

1043 

261 

1451 

FNS 

425 

197 

54 

262 

STEPS 

2805 

897 

3202 

1390 

JACS 

ir^ 

15783 

6764 

17954 

7923 

FNS 

2834 

1443 

3245 

1423 

STEPS 

2683 

1383 

4655 

6645 

JACS 

lOr* 

15160 

10362 

26317 

37297 

FNS 

2726 

•  2269 

4789 

6760 

STEPS 

Table  4.3b.  The  number  of  Jacobian  evaluations  (JACS),  function  evalua- 
tkMS  (FNS),  and  time  steps  (STCPS)  needed  to  solve  ptoblmi  FS  from  the 
stiff  test  set  with  DASSL  and  finite-diffioence  Jacobians  using  the  standaid. 
STAB,  and  PI  controllers. 


